Acces parameter recovery results for PH with C model

1 Setup

rm(list=ls())
libs<-c("rstan", "gdata", "bayesplot", "stringr", "dplyr", "ggplot2", "hBayesDM")
sapply(libs, require, character.only=TRUE)

Loading required package: rstan

Loading required package: StanHeaders

rstan version 2.26.22 (Stan version 2.26.1)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

Loading required package: gdata

Attaching package: 'gdata'

The following object is masked from 'package:stats':

    nobs

The following object is masked from 'package:utils':

    object.size

The following object is masked from 'package:base':

    startsWith

Loading required package: bayesplot

This is bayesplot version 1.10.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

Loading required package: stringr

Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:gdata':

    combine, first, last

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Loading required package: ggplot2

Loading required package: hBayesDM

Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called 'hBayesDM'

    rstan     gdata bayesplot   stringr     dplyr   ggplot2  hBayesDM 
     TRUE      TRUE      TRUE      TRUE      TRUE      TRUE     FALSE 

datapath <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/ILT_Stan_Modeling'
out_path <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/ILT_Stan_Modeling/Output'
behavpath <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/ILT_DATA'

# load files containing true parameters used as input for simulation
orig_file <- 'fit_n_142_2024-01-29_bandit2arm_delta_PH_withC_estimation1_delta0.9_stepsize0.1.rds'
orig_fit <- readRDS(file.path(out_path, orig_file)) # Stan model output

# load file containing true parameters computed as transformed parameters during simulation
sim_file <- 'sim_2024-01-31_bandit2arm_delta_PH_withC_sim.rds'
sim_fit <- readRDS(file.path(out_path, 'Parameter_Recovery', sim_file))

# get file names for simulated data results
recovery_file <- 'recovery_2024-02-01_bandit2arm_delta_PH_withC.rds'
recovery_fit <- readRDS(file.path(out_path, 'Parameter_Recovery', recovery_file)) # Stan model output
  
color_scheme_set("mix-blue-pink")

2 Posterior predictive checks & correlations

# Load true parameters

## extract posterior means for all parameters to use them as input for simulation

### posterior means of parameters as input for simulation
true_mu_pr <- as.vector(summary(orig_fit, pars="mu_pr")$summary[, c("mean")]) 
true_sigma <- as.vector(summary(orig_fit, pars="sigma")$summary[, c("mean")]) 

true_A_pr <- as.vector(summary(orig_fit, pars="A_pr")$summary[, c("mean")]) 
true_tau_pr <- as.vector(summary(orig_fit, pars="tau_pr")$summary[, c("mean")]) 
true_gamma_pr <- as.vector(summary(orig_fit, pars="gamma_pr")$summary[, c("mean")]) 
true_C_pr <- as.vector(summary(orig_fit, pars="C_pr")$summary[, c("mean")]) 

### transformed parameters, output of simulation
true_sim_pars <- as.vector(extract(sim_fit))
true_A <- as.vector(true_sim_pars$A)
true_tau <- as.vector(true_sim_pars$tau)
true_gamma <- as.vector(true_sim_pars$gamma)
true_C <- as.vector(true_sim_pars$C)

true_mu_A <- as.vector(true_sim_pars$mu_A)
true_mu_tau <- as.vector(true_sim_pars$mu_tau)
true_mu_gamma <- as.vector(true_sim_pars$mu_gamma)
true_mu_C <- as.vector(true_sim_pars$mu_C)
  
## extract parameter values based on simulated data
recovered_mu_pr <- as.matrix(recovery_fit, pars = "mu_pr")
recovered_sigma <- as.matrix(recovery_fit, pars = "sigma")

recovered_A_pr <- as.matrix(recovery_fit, pars = "A_pr")
recovered_tau_pr <- as.matrix(recovery_fit, pars = "tau_pr")
recovered_gamma_pr <- as.matrix(recovery_fit, pars = "gamma_pr")
recovered_C_pr <- as.matrix(recovery_fit, pars = "C_pr")

recovered_A <- as.matrix(recovery_fit, pars = "A")
recovered_tau <- as.matrix(recovery_fit, pars = "tau")
recovered_gamma <- as.matrix(recovery_fit, pars = "gamma")
recovered_C <- as.matrix(recovery_fit, pars = "C")

recovered_mu_A <- as.matrix(recovery_fit, pars = "mu_A")
recovered_mu_tau <- as.matrix(recovery_fit, pars = "mu_tau")
recovered_mu_gamma <- as.matrix(recovery_fit, pars = "mu_gamma")
recovered_mu_C <- as.matrix(recovery_fit, pars = "mu_C")

## Compare true and recovered parameters
mcmc_recover_intervals(recovered_mu_pr, true_mu_pr, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_mu_pr, true_mu_pr)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_mu_A, true_mu_A, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_mu_A, true_mu_A)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_mu_tau, true_mu_tau, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_mu_tau, true_mu_tau)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_mu_gamma, true_mu_gamma, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_mu_gamma, true_mu_gamma)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_mu_C, true_mu_C, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_mu_C, true_mu_C)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_sigma, true_sigma, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_sigma, true_sigma)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_A_pr, true_A_pr, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_A_pr, true_A_pr)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_tau_pr, true_tau_pr, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_tau_pr, true_tau_pr)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_gamma_pr, true_gamma_pr, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_gamma_pr, true_gamma_pr)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_C_pr, true_C_pr, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_C_pr, true_C_pr)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_A, true_A, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_A, true_A)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_tau, true_tau, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_tau, true_tau)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_gamma, true_gamma, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_gamma, true_gamma)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_recover_intervals(recovered_C, true_C, prob = 0.5, prob_outer = 0.95)

mcmc_recover_hist(recovered_C, true_C)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.